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Abstract. The Antarctic bedrock is evolving as the solid 
Earth responds to the past and ongoing evolution of the ice 
sheet. A recently improved ice loading history suggests that 
the Antarctic Ice Sheet (AIS) has generally been losing its 
mass since the Last Glacial Maximum. In a sustained warm- 
ing climate, the AIS is predicted to retreat at a greater pace, 
primarily via melting beneath the ice shelves. We employ the 
glacial isostatic adjustment (GIA) capability of the Ice Sheet 
System Model (ISSM) to combine these past and future ice 
loadings and provide the new solid Earth computations for 
the AIS. We find that past loading is relatively less important 
than future loading for the evolution of the future bed to- 
pography. Our computations predict that the West Antarctic 
Ice Sheet (WAIS) may uplift by a few meters and a few tens 
of meters at years AD 2100 and 2500, respectively, and that 
the East Antarctic Ice Sheet is likely to remain unchanged 
or subside minimally except around the Amery Ice Shelf. 
The Amundsen Sea Sector in particular is predicted to rise 
at the greatest rate; one hundred years of ice evolution in this 
region, for example, predicts that the coastline of Pine Is- 
land Bay will approach roughly 45mmyr -1 in viscoelastic 
vertical motion. Of particular importance, we systematically 
demonstrate that the effect of a pervasive and large GIA up- 
lift in the WAIS is generally associated with the flattening 
of reverse bed slope, reduction of local sea depth, and thus 
the extension of grounding line (GL) towards the continen- 
tal shelf. Using the 3-D higher-order ice flow capability of 
ISSM, such a migration of GL is shown to inhibit the ice 
flow. This negative feedback between the ice sheet and the 


solid Earth may promote stability in marine portions of the 
ice sheet in the future. 


1 Introduction 

Projecting the evolution of the Antarctic Ice Sheet (AIS) 
into the next few centuries relies on simulating a complex 
and non-linear coupled Earth system. A recent survey of ex- 
perts by Bamber and Aspinall (2013) reveals that projections 
for AIS contribution to the rate of sea level rise at the year 
AD 2100 are generally rather moderate (~ 1.7 mmyr -1 ) and 
that the upper end of the spectrum of projections would be 
about 7 times this value, mainly owing to intensification of 
dynamics of the West Antarctic Ice Sheet (WAIS). However, 
projections beyond AD 2100 are much more uncertain (Bind- 

schadler et al., 2013) and are mainly limited by the poor 
knowledge of the physics involved in the grounding line (GL) 
migration and ice shelf melting and calving (e.g., Vaughan 
and Arthern, 2007; Walker et al., 2013). Furthermore, there 
is strong evidence that over the past four million years, dur- 
ing times of increased global atmospheric temperatures by 
2-3 °C, the marine WAIS collapsed, and possibly some por- 
tions of the larger East Antarctic Ice Sheet (EAIS) as well 
(Naish et al., 2009; Raymo et al., 2011; Cook et al., 2013). 

If we are to improve our abilities to assess the risk of 
the catastrophic consequences of partial collapse of AIS 
marine-based ice currently locked to the Antarctic continent, 
steps must be taken to assess fully the role of solid Earth 
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deformation over tens to thousands of years, during which 
time gravitational viscoelastic flow of the underlying man- 
tle may act to change the stability of the AIS. Past assess- 
ments have been that isostatic uplift following the ice sheet 
retreat promoted stability of the WAIS near the end of the 
last deglaciation during the mid-Holocene (Thomas, 1976). 
More recently, the promotion of stability by glacial isostatic 
adjustment (GIA) has been shown with increasingly sophisti- 
cated computations (Gomez et ah, 2010, 2013). It is also now 
recognized that past AIS recession and ice flow direction are 
plausibly explained by strong interaction of ice loading with 
the solid Earth (Siegert et al., 2013). It is this ice-sheet/solid- 
Earth (IS/SE) interaction that we now explore in this paper 
for the AIS as a whole, using the vertical surface motions 
of bedrock GIA (Ivins and James, 1999) and the 3-D ther- 
momechanical ice flow (Pattyn, 2003) capability of the Ice 
Sheet System Model (ISSM) (Larour et al., 2012b). 

1.1 IS/SE feedback: GL migration 

Perhaps the most important IS/SE feedback is associated 
with the GL migration (e.g., Lingle and Clark, 1985). For 
equilibrium sea level, any change in vertical bedrock eleva- 
tion also changes the local sea depth. This perturbs buoyant 
forces in the regional ocean water and may promote the mi- 
gration of the GL. Uplift of the seabed occurs, for example, 
in response to thinning of the inland ice sheet, and it causes 
local sea depth to decrease. If the GL is in initial equilibrium 
but sea depth decreases due to bed uplift, the GL tends to 
advance towards the continental shelf (e.g., Weertman, 1974; 
Schoof, 2007). The conceptual illustration of this negative 
feedback is shown in Fig. la. Lingle and Clark (1985) ex- 
plored the effect of GIA-related seabed uplift on GL mi- 
gration for Ice Stream E, now known as the MacAyeal Ice 
Stream, in the WAIS. The modeled GIA uplift, caused by 
thinning of the ice stream catchment area, delays the onset of 
GL retreat, thus reducing the rate of retreat during Holocene 
sea level rise. Notably, it was argued that the regional ad- 
vance of Ross Sea GL may have occurred over the past three 
millennia. The gravitational attraction effects on local sea 
level developed later by Gomez et al. (2010) act to amplify 
this negative feedback during ice sheet retreat, since the di- 
minished mass behind the GL has less mutual gravitational 
attraction with adjacent sea water, thus causing local sea level 
to drop. 

The pace and magnitude of GL migration are also dictated 
by the bedrock slope. If sea depth decreases (due to bed up- 
lift and sea level drop associated with the thinning of inland 
ice) at the initially stable GL on a reverse bed slope, for ex- 
ample, the GL tends to advance further on a relatively flat 
bed. Therefore the GLs associated with the Ross and Ronne 
ice shelves in the WAIS having relatively flat reverse bed 
slopes are sensitive to this effect (Conway et al., 1999). Cap- 
turing the dynamics of such GL sensitivity demands proper 
understanding of interactions between the ice, the ocean and 


the solid Earth, and is indeed the key to successful model- 
ing of the WAIS (e.g., Vaughan and Arthern, 2007; Katz and 
Worster, 2010). 

1.2 Additional IS/SE feedbacks 

There are other important feedback mechanisms that the 
solid Earth deformation provides to ice sheets. The GIA 
uplift can be important in providing basal resistance to ice 
flow and buttressing the ice sheet by raising bedrock pinning 
points (e.g., Favier et al., 2012; Siegert et al., 2013). This 
concept is captured in Fig. lb. The GIA-induced changes in 
surface elevation and regional slope of the ice sheet may af- 
fect the gravitational driving stress, as well as some processes 
at the ice-atmosphere interface (e.g., surface mass balance). 
These perturb the momentum balance and affect the englacial 
velocity field (e.g., Cuffey and Paterson, 2010; Winkelmann 
et al., 2012), which in turn may impact the ice thermody- 
namics via changes in strain heating. Spatially varying bed 
uplift also affects the hydraulic potential field and hence the 
subglacial hydrology and the sliding rate of the ice sheet. 
There is additionally the complication of bulge migration, 
a broad-scale phenomenon involving bending of the elastic 
lithosphere and mantle lateral flow. Due to the lateral motion 
of this topographic bulge, local crustal motions (and slopes) 
may change sign during GIA (e.g., Fjeldskaar, 1994). 

These mechanisms are extremely difficult to isolate and 
quantify, and it is therefore not obvious whether (and un- 
der what circumstances) each of these acts to accelerate or 
inhibit the ice flow. As long as the thermomechanical ice 
sheet model and other companion models (e.g., surface mass 
balance model and the hydrological model) are dynamically 
coupled to a comprehensive solid Earth model, however, 
most of these feedbacks are intrinsically taken into account. 

1.3 New solid Earth computations 

For large timescale (millennia) simulations, most of the ice 
sheet models (e.g., Pollard and DeConto, 2009; Hughes et al., 
2011) capture the solid Earth physics of varying degrees 
of complexity (cf. Le Meur and Huybrechts, 1996). With 
the exception of the recent work of Gomez et al. (2013), 
none of these studies provides the explicit assessment of 
effects of GIA uplift on several aspects of ice sheet dy- 
namics (e.g., GL migration and gravitational driving stress). 
Le Meur and Huybrechts (2001) explicitly pointed out the 
need for the more complete coupling that could be found in 
the wavelength-dependent relaxation spectra of viscoelastic 
solid Earth models. In this paper, we quantify two distinct 
IS/SE feedback mechanisms applied to the AIS on centennial 
timescales using multiple wavelength-dependent decay spec- 
tra. Assuming the equilibrium sea level, we first evaluate how 
the future bed uplift (governed by the past and future evolu- 
tion of AIS; cf. Sect. 2) controls the GL migration. Second, 
we assess the role of bed uplift in the gravitational driving 
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Figure 1 . Schematic of IS/SE feedbacks, (a) GL migration controlled by the local sea depth. For equilibrium sea level, the GIA uplift 
due to thinning of the inland ice promotes the GL migration towards the continental shelf. The hydrostatic equilibrium requires that s\/b\ = 
s l/^2 — — (Pw~ P)/P, where S[ and bi are the surface and bedrock elevations (subscript i = 1,2 represents the initial and final configurations, 
respectively), and p w and p are the water and ice densities, (b) Pinning point, raised by GIA uplift due to thinning of the inland ice, provides 
basal resistance and buttressing to the ice sheet. Red arrows depict the velocity profiles. Dashed lines in the final configurations (second 
column) represent the initial geometry. Note that the destabilizing effects of GIA during the periods of inland ice thickening that causes the 
bedrock to subside can be conceptualized by reverting the direction of mid-arrows. 


stress. The overall influence of solid Earth deformation (via 
changes in GL and driving stress) on the ice surface veloci- 
ties is also quantified. These assessments, based on reliable 
models and data, help us to understand whether the GIA ef- 
fects are significant in controlling the future evolution of AIS 
on centennial timescales. 

The paper layout is as follows. In Sect. 2, we introduce the 
solid Earth model employed in this research, discuss how the 
change in ice thickness is translated into the ice load height, 
provide a detailed account of model tuning procedure, and 
describe all of the required model input data. In Sect. 3, we 
present modeling results of the future bedrock topography, 
and assess the relative importance of the past and future ice 
load changes. In Sect. 4, we quantify the influence of pre- 
dicted change in Antarctic bed topography on several aspects 
of ice sheet dynamics. Finally, in Sect. 5, we summarize key 
conclusions of broad interest. 


2 Model and data 

Ivins et al. (2013) presented a much improved ice loading 
history for the AIS since the Last Glacial Maximum (LGM). 
As for the future, the recently concluded SeaRISE (Sea-level 
Response to Ice Sheet Evolution) project (Bindschadler et al., 
2013; Nowicki et al., 2013) provided quantitative projections 
of the evolution of AIS under ongoing climate warming. We 
employ the new GIA capability (Ivins and James, 1999) of 


ISSM (Larour et al., 2012b), hereinafter referred to as the 
ISSM/GIA model, to combine these data of past and future 
ice loadings and calculate the change in bedrock topography 
over the same timescale projections as in the SeaRISE stud- 
ies. Using appropriate analytical and numerical models, we 
then evaluate the stabilizing or destabilizing influence of pre- 
dicted changes in bed topography on the ice sheet dynamics. 

While the SeaRISE experiments employed state-of-the-art 
numerical treatments of ice flow, it should be noted that the 
majority of these models were not coupled to a comprehen- 
sive solid Earth model. Furthermore, the SeaRISE experi- 
ments do not capture the paleo-evolution of the AIS since 
the LGM, thus limiting the possibility for the participating 
ice sheet models (with GIA capability) to perform similar 
analysis presented in this study. 

2.1 The ISSM/GIA model 

ISSM is a continental-scale, high-resolution, multi-model 
simulation code developed for understanding the dynamic 
behavior of large ice sheets (Larour et al., 2012b). This 
open source finite-element software is capable of simulat- 
ing ice-flow mechanics of varying degrees of complexity 
(Seroussi et al., 2012), performing sensitivity analysis (e.g., 
Larour et al., 2012a), inverting unknown field parameters 
(e.g., Morlighem et al., 2010), and assessing mechanics of 
rift propagation and eventual collapse of ice shelves (e.g., 
Borstad et al., 2012). The semi-analytical GIA solution of 
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Ivins and James (1999) is one of several new features being 
actively implemented in ISSM. Here we briefly summarize 
the key elements of the model. 

We assume that the ice sheet rests on top of the solid non- 
rigid Earth, which is considered to be a simple two-layered 
incompressible continuum with an upper elastic lithosphere 
floating on the viscoelastic (Maxwell material) mantle half- 
space. The theory governing the deformation of pre-stressed 
solid Earth subject to a normal surface traction of ice sheet 
relies upon the fundamental equations of motion and is dis- 
cussed elsewhere (e.g., McConnell, 1965; Wolf, 1985; Ivins 
and James, 1999). For axisymmetric loading problems, it is 
possible to obtain the semi-analytical solution of vertical dis- 
placement at the lithosphere surface (i.e., both the ice-bed 
and ocean-bed interfaces). This is the essence of the solu- 
tion for GIA assessment. For the equilibrium sea level hy- 
pothesis, however, the GIA solutions perturb the ice sheet 
only within the area of grounded ice. For the AIS that is sur- 
rounded on most of its periphery by floating ice, the extent of 
the grounded ice may evolve, as we assume in this study, ac- 
cording to the hydrostatic balance between the ice shelf and 
ocean water. 

Given the appropriate ice loading history and choice of the 
model/material parameters (cf. Sect. 2.3), the GIA solution 
depends on the size of the ice load itself and the radial dis- 
tance of the evaluation point from the load center. Assumed 
axisymmetry implies that the shape of the ice load be essen- 
tially cylindrical (e.g., Ivins and James, 1999). In the Carte- 
sian framework of the ISSM/GIA model, we treat the size of 
the ice load as the property of a mesh element and compute 
the GIA solution at each node of the element. Individual 2- 
D (xy plane) mesh elements are defined as the equivalence 
of a footprint (i.e., projection onto the vy plane) of cylindri- 
cal disc loads, ensuring that the corresponding element and 
disc both share the same origin and plan- form area (cf. inset 
of Supplement Fig. SI a). The height of the ice load is then 
assigned to each element such that the average normal trac- 
tional force on the corresponding area of ice-bedrock contact 
is conserved. At each node within the domain, the final GIA 
solutions are computed by integrating the solutions due to in- 
dividual disc loads, defined as the property of mesh elements. 

The ISSM/GIA model is tested against the benchmark ex- 
periments (Wolf, 1985; Ivins and James, 1999) and found to 
be sensitive to the mesh resolution. For reasonably fine reso- 
lution (element size typically on the same order of magnitude 
as ice thickness, or two orders of magnitude smaller than the 
characteristic size of an ice sheet), however, the model per- 
forms well within the acceptable accuracy. Sample results are 
provided in Supplement Fig. SI. 

2.2 Differential ice height 

Prior to describing the ISSM/GIA model applied to the AIS, 
we briefly discuss how the change in Antarctic ice thickness 
is translated into the height of ice load. In this study, we as- 


sume that the sea level remains in its present-day state. We 
define the height of ice load at any time, t , as the differential 
ice height (DIH), A h(x,y,t), with reference to the present- 
day configuration of the AIS (James and Ivins, 1998). 

In the regions where ice is grounded both presently and at 
the given time t, DIH is simply the change in ice thickness 
over the corresponding time period. Similarly, in the regions 
where ice is floating both at present and at time t, DIH is 
zero as we assume that ice shelves are freely floating over the 
ocean, respecting the hydrostatic equilibrium. Defining DIH 
is a bit complicated in the areas over which the GF migrates 
over the course of time. If ice thickness at time t is smaller 
than the present-day value (i.e., ice floats at time t , but is now 
grounded), the DIH is defined as follows: 

A h(x, y, t) = — b(x, y, t ) — h(x, y, 0), (1) 

P 

where p w is the ocean water density, p is the ice density, 
b(x , y, t) < 0 is the isostatically corrected bedrock elevation 
(with respect to the present-day sea level) of the marine por- 
tions at time t (James and Ivins, 1998), and h(x, y, 0) is the 
present-day ice thickness. DIH in such a case would be neg- 
ative. Similarly, in the areas where the ice sheet holds thicker 
ice at time t than at present (i.e., ice now floating is grounded 
at time t), DIH is defined as follows: 

A h(x, y, t) = h{x, y, t) — — b(x, y, 0), (2) 

P 

where h(x, y, t) is the ice thickness at time t and b(x,y, 0) < 
0 is the present-day bedrock elevation of the marine portions 
of the ice sheet. DIH in such a case would be positive. 

As a general convention in this study, we use t < 0 to de- 
note the past, i.e., before present, and t > 0 for the future 
unless stated otherwise. 

2.3 Model tuning 

We apply the ISSM/GIA model to the AIS. We mesh the foot- 
print of present-day AIS (Bamber et al., 2009) into triangular 
elements. In order to capture the potentially interesting fea- 
tures, the domain is discretized to consist of a high-resolution 
mesh around the coast (typical element size of 1 0 km) than 
in the interior of the ice sheet (element size of 25 km). This 
unstructured mesh captures the model inputs (e.g., past or 
future ice loads) in sufficient detail and provides a reason- 
able compromise between solution accuracy and computa- 
tional efficiency. Doubling the mesh density, for example, 
improves the GIA solution (under present-day ice loading) 
only slightly (< 0.1 %), but it increases the high computa- 
tional cost by one order of magnitude. 

A comprehensive list of modern global positioning system 
(GPS) measurements is presented by Thomas et al. (2011). 
However, we tune our model by testing it against 18 high- 
precision data. Following Ivins et al. (2013), we first elimi- 
nate records from the Antarctic Peninsula north of 72° S due 
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to the associated difficulty of dealing with large elastic and 
transitional viscoelastic signals present there. We then aver- 
age the values from stations located within 100 km of one an- 
other and eliminate some stations with reported errors greater 
than the signal amplitude. 

In order to make reasonable predictions of the present-day 
GIA uplift, the slow response of highly viscous solid Earth 
demands that the evolution of AIS during the past several 
thousand years be considered in the ISSM/GIA model. There 
are a few reliable GIA ice loading histories for Antarctica 
(e.g., Peltier, 2004; Ivins and James, 2005; Whitehouse et al., 
2012). These generally describe the timing and magnitude of 
deglaciation since the LGM based on geological and ice core 
data. In this study, we employ a much improved loading his- 
tory discussed in Ivins et al. (2013). By upgrading the load- 
ing history of Ivins and James (2005) with recently available 
geochronological constraint data, the later model was able to 
provide a more accurate GIA correction to GRACE (Gravity 
Recovery and Climate Experiment) data measured over the 
period AD 2003-2012. 

Based on Ivins et al. (2013), we have isostatically cor- 
rected DIHs available for 1 1 time stamps in the past (at — 1 , 
-2.2, -3.2, -6.8, -7.6, -11.5, -15, -17, -19, -21 and 
— 102kyr; see Supplement Fig. S2). Note that t = — 21 kyr 
roughly corresponds to the LGM of the AIS, while —15 kyr 
marks the onset of deglaciation of the WAIS (e.g., Clark 
et al., 2009). For t = — 1 kyr, we consider zero DIH that could 
be constrained using the recently available surface mass bal- 
ance data (e.g., Verfaillie et al., 2012; Favier et al., 2013). 
However, this process is not straightforward because the 
magnitude and spatial distribution of ice flux during the pe- 
riods of inferred mass balance are vastly unknown. Further- 
more, we consider t = —102 kyr to mark the initial configu- 
ration for AIS as being identical to the present-day configura- 
tion. This implicitly assumes that the DIHs before the LGM 
have a minimal impact on the current and future response of 
the solid Earth. (We demonstrate in Sect. 3.2 that this is in- 
deed a valid assumption.) Note also that the ice loading on 
the ISSM/GIA model is assumed to vary in a piece-wise lin- 
ear fashion between the adjacent time stamps. 

The model and material parameters considered in this 
study approximate the preliminary reference Earth model 
(Dziewonski and Anderson, 1981) and are taken from Ivins 
et al. (2013) unless otherwise specified. For several val- 
ues for lithosphere thickness, Ivins et al. (2013) performed 
a parameter-space study in their two-layer mantle model. 
Noting that 65 and 115 km may represent, respectively, the 
average thickness of the lithosphere for the WAIS and EAIS, 
Ivins et al. (2013) provided appropriate combinations of the 
upper and lower mantle viscosity. Because the past (see Sup- 
plement Fig. S2) and future DIHs (cf. Sect. 2.4 and Sup- 
plement Fig. S3) indicate that the majority of changes oc- 
cur in the WAIS, we consider a 65 km thick lithosphere in 
our model. We cannot pick the corresponding mantle viscos- 
ity from Ivins et al. (2013), as our model does not have two 


mantle layers. We therefore consider the mantle viscosity as 
a tuning parameter, such that the difference in the mean be- 
tween the measured modern GPS data (Thomas et al., 2011) 
and modeled current GIA uplift at 18 data points is min- 
imized (Fig. 2b). The optimized solutions for current up- 
lift rate are shown in Fig. 2a. Key characteristic features 
of our predictions include greater uplift rate around the Mt. 
Ellsworth territory and a mild rate of bed subsidence in 
the interior of EAIS. Such spatial patterns of uplift rate es- 
sentially reflect signatures of the employed ice loading his- 
tory (cf. Supplement Fig. S2). Note that the optimal pre- 
dictions of uplift rate (Fig. 2a) correspond to a mantle vis- 
cosity of 7 x 10 20 Pas. As expected, this magnitude falls in 
between the upper (2 x 10 20 Pas) and lower mantle viscos- 
ity (1 .5 x 10 21 Pa s) recommended for the chosen lithosphere 
thickness (Ivins et al., 2013). While the architecture of the 
ISSM/GIA model can capture high-resolution spatial vari- 
ability in solid Earth material parameters, we do not experi- 
ment with lateral inhomogeneities in this study. 

2.4 Future ice loading 

The AIS mass change may become more dynamic in the fu- 
ture due to ice shelf melting (e.g., Pritchard et al., 2012; De- 
poorter et al., 2013; Rignot et al., 2013). The SeaRISE partic- 
ipating ice sheet models, primarily driven by melt-dominated 
forcing, quantified the future evolution of AIS under the 
so-called “R8 scenario”, which is the proxy of representa- 
tive concentration pathway emission scenario 8.5 (RCP 8.5) 
(Bindschadler et al., 2013; Nowicki et al., 2013). The RCP 

8.5 scenario represents an ongoing rise in emissions through- 
out the century, reaching 8.5 Wm -2 at AD 2100 (e.g., van 
Vuuren et al., 2011). The radiative forcing associated with 
the RCP 8.5 scenario is loosely related to the R8 scenario via 
all three components of the SeaRISE model forcing, namely 
the surface climate, basal sliding and ice shelf melting. As 
these forcings are the ones that govern the future evolution 
of AIS, it is relevant in the present context to provide a brief 
account of them. 

Firstly, the SeaRISE surface climate forcing follows a 

1.5 x A1B scenario (IPCC-AR4, 2007) until AD 2200. (The 
A IB scenario generally describes a future world of rapid 
growth in economy and population that peaks in mid-century, 
and technologies that rely equally on both fossil and non- 
fossil sources of energy.) A mild increase in surface tempera- 
ture, a total of 0.5 °C, is assumed during the period AD 2200- 
2500. Secondly, no sliding amplification is considered until 
AD 2100 assuming that the Antarctic surface temperature 
will remain below zero, thus ignoring the potential for sur- 
face melt-induced basal sliding prior to this time. Thereafter, 
sliding increases linearly at a rate of 20 % of its original value 
per century, but only in coastal regions. Inland, the sliding 
amplification factor decreases linearly as a function of sur- 
face elevation such that no sliding enhancement is applied 
above 1200ma.s.l. Thirdly, ice shelf melting is assumed to 
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Figure 2. Model tuning and predictions of current uplift rate, (a) Modeled GIA uplift rate at the present day. Calculations are made by 
forcing the ISSM/GIA model by ice loading history over the past 21 kyr (Supplement Fig. S2) (Ivins et al., 2013). Black circles locate the 
position where model results are within the 1-cr uncertainty range of GPS measurements (Fig. 2b). Red circles represent cases where the 
model overestimates the measurements, whereas blue circles indicate locations where the model underperforms. Big circles are to denote the 
absolute misfits that are >0.75 mmyr -1 . (b) Validation of the model against 18 high-precision GPS uplift data (Thomas et al., 2011). Error 
bars depict 1-a uncertainties associated with the GPS measurement. 




increase linearly from its present-day value to 60myr _1 at 
AD 2200. Additional lOrnyr -1 melt extends linearly over 
the next 300 yr. Changes in basal melting conditions are only 
applied to the Amundsen Sea Sector (90-120° W) and the 
Amery Ice Shelf (60-75° E), not to the Weddell and Ross 
seas. Such a restriction of ice shelf melting considered in the 
SeaRISE experiment seems reasonable with reference to cur- 
rent observations (Depoorter et al., 2013; Rignot et al., 2013) 
that reveal that ice shelves around the Amundsen Sea are the 
most susceptible to melting. Furthermore, spatial distribution 
of ocean temperature anomalies (e.g., Pritchard et al., 2012) 
also supports this hypothesis. 

Under the R8 scenario, a total of four ice sheet models 
simulates the future evolution of the AIS through AD 2500. 
These models are the Anisotropic Ice Flow Model (AIF) 
(Wang et al., 2012), the Penn State Ice Sheet Model 
(PennState) (Pollard and DeConto, 2009), the Potsdam Par- 
allel Ice Sheet Model (PISM-PIK) (Winkelmann et al., 
2011), and the Simulation Code for Polythermal Ice Sheets 
(SiCoPolIS) (Greve, 1997). These ice sheet models em- 
ploy different assumptions and methods for solving the full 
physics involved in simulating ice flow (e.g., shallow-ice 
vs. first-order flow mechanics; different treatments for basal 
sliding, subglacial hydrology, and GL migration), they em- 
ploy different numerics (e.g., different spatial and temporal 
resolutions), they have unique techniques for dealing with 
the prescribed model forcing, and so forth (cf. Bindschadler 
et al., 2013; Nowicki et al., 2013). Consequently, each ice 
sheet model produces a unique evolution of the AIS. We ex- 
tract the ice thickness from each model prediction for five 
time stamps (at t = 100, 200, 300, 400 and 500 yr) and com- 
pute DIHs using Eqs. (1) and (2) as appropriate. For simplic- 
ity, we assume b(x,y,t)~b(x,y,0) in Eq. (1) (i.e., future 
bedrock elevation is not isostatically corrected while com- 


puting future DIHs). We anticipate that the associated er- 
rors translated into the GIA solutions would be minimal. Ex- 
amples of future DIHs are provided in Supplement Fig. S3. 
Again, ice loads in the ISSM/GIA model are assumed to vary 
linearly between the adjacent time stamps. 

3 Future bed topography 

In order to predict the future bed topography for AIS, the 
calibrated ISSM/GIA model (cf. Sect. 2.3) is forced by an 
appropriate sequence of ice load changes into the future. 
As we have four independent sets of future ice loading (cf. 
Sect. 2.4), we may compute four unique GIA solutions at 
any evaluation time in the future. Based on these solutions, 
here we present the estimates of future bed uplift, isolate the 
role of past and future ice loading, and also evaluate how 
predicted change in bed topography alters the bedrock slope. 
Note that in this section we deviate from the general time 
convention and present our results (predictions) in “AD”. 

3.1 Vertical bed displacement 

The GIA solutions for individual future ice loading, com- 
bined with the consideration of a lone spin-up loading his- 
tory, are computed at AD 2100 and AD 2500 and shown in 
Supplement Figs. S4 and S5, respectively. Although these so- 
lutions differ from each other in both magnitude and their 
spatial distribution, some common features are noteworthy: 
(i) all models predict minor subsidence in a few places, par- 
ticularly along the Wilkes Land coast; (ii) the topography 
of the interior of the EAIS is likely to remain unchanged; 
and (iii) a pervasive and large uplift is predicted in the WAIS 
(except for the AIF simulations) and around the Amery Ice 
Shelf (except for the PennState and PISM-PIK simulations). 
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We can, of course, offer no one GIA solution as being 
more reliable than any other. We therefore calculate the av- 
erage of model solutions (hereafter termed “model-average 
solutions”) (Fig. 3a and b) and consider these as reason- 
able estimates of the Antarctic bed uplift. It is possible that 
our ensemble approach for these predictions is insufficient. 
Nonetheless, we assert that they provide the correct order- 
of-magnitude estimates and the likely spatial patterns of the 
future bed uplift, which are sufficient to evaluate some of the 
IS/SE interactions outlined in Sect. 1.3. 

By AD 2100, the Amundsen Sea Sector and the Amery 
Ice Shelf may rise by about four and three meters, respec- 
tively (Fig. 3a). The rest of the WAIS is likely to rise by 
up to two meters. The interior of the EAIS is predicted to 
remain unchanged. The Adelie and Wilkes lands, where all 
ice sheet models predict large snow accumulation (Nowicki 
et al., 2013), are likely to subside by about less than one me- 
ter. It should be noted that, for the chosen climate change 
scenario, all ice sheet models but SiCoPolIS predict moder- 
ate snow accumulation in the Queen Maud, Ross and Wed- 
dell basins, and minimal accumulation in the Amundsen and 
Amery basins (cf. Nowicki et al., 2013). Roughly similar spa- 
tial patterns of GIA uplift are predicted for AD 2500 as well 
(Fig. 3b). In this case, the bed may uplift by about 25 m in the 
Amundsen Sector, and by about 10-15 m in the rest of WAIS 
and the Amery Ice Shelf. The interior of EAIS may remain 
mostly unperturbed. Bed subsidence of about four meters is 
likely to occur along the Adelie and Wilkes Land coasts, as 
well as along the coast north of the Amery Ice Shelf. 

We also find similar spatial patterns for the model-average 
bed uplift rates (cf. Supplement Fig. S6). At AD 2100, 
the Amundsen Sea Sector and Wilkes Land are predicted 
to rise and subside at the highest rates of about 45 and 
— 5 mmyr -1 , respectively. The interior of Marie Byrd Land 
is predicted to rise at the highest rate of 75 mmyr -1 at 
AD 2500; The Amundsen Sea Sector also rises at a large 
rate of about 40 mmyr -1 . The greatest rate of subsidence (of 
about —15 mmyr -1 ) is predicted along the east coast except 
around the Amery Ice Shelf. 

At both evaluation times, as noted earlier, we obtain differ- 
ent GIA solutions associated with individual future ice load- 
ing. Here we briefly outline how well these predictions for 
the Antarctic bed uplift match one another. The standard de- 
viation shown in Fig. 3c illustrates that the model predictions 
are generally in good agreement with each other at AD 2100, 
except around the Amundsen Sector and the Amery Ice 
Shelf. As depicted in Fig. 3d, similar agreement is found 
amongst the model predictions at AD 2500. Large deviations 
are predicted once again around the Amundsen Sea Sector. 
Moderate deviations can be seen along the east coast includ- 
ing the Amery Ice Shelf. Such deviations amongst model pre- 
dictions for the Antarctic bed topography predicted both at 
AD 2100 and AD 2500 can generally be attributed to the 
limiting values of GIA solutions predicted by the PennState 
and AIF models (cf. Supplement Figs. S4 and S5). In both 


years, PennState predicts large uplift in the Amundsen Sec- 
tor and the least uplift (in fact, subsidence) around the Amery 
Ice Shelf; the opposite is true in the case of AIF predictions. 

Our predictions might slightly overestimate the GIA so- 
lutions in the EAIS, as the modeled lithosphere thickness 
(i.e., 65 km) is much thinner than the more common value 
(115 km, as reported in Sect. 2.3). Nevertheless, the follow- 
ing general findings should remain unaltered. With refer- 
ence to the WAIS GIA solutions, (i) Amery Ice Shelf may 
rise moderately; (ii) the interior of the EAIS may remain 
mostly unperturbed; and (iii) minor subsidence may occur 
along the coastal EAIS. It also needs to be mentioned here 
that we do not model the subglacial erosion in this study, al- 
though quite rapid evacuation of soft sediments is now occur- 
ring at the bed of the Pine Island Glacier. Erosion has been 
roughly estimated to cause the topography to lower at a rate 
of 0.6 zb 3 myr -1 as ice flow exceeding 1 kmyr -1 erodes ma- 
terial in deep longitudinal fjord valleys (Smith et al., 2012). 
This erosive action typically takes place in 20 km wide val- 
leys of approximately 200 km length. While an important 
consideration in ice sheet modeling (Kessler et al., 2008), 
GIA topographic responses occur over much broader length 
scales exceeding the areal dimensions of fast erosion by 
nearly two orders of magnitude, thus having an impact on 
the evolution of the ice sheet on the scale of the drainage 
basin itself. 

3.2 Role of past and future loading 

Our predictions of bed uplift reflect the combined effects 
of long-term viscous creep of solid Earth driven by the ice 
loading history and its short-term (centennial timescale) vis- 
coelastic response to the future ice load change. It might be 
useful to isolate the contribution of past and future ice load- 
ing to the evolution of future bed topography. First, we let 
the calibrated ISSM/GIA model (driven by the past load- 
ing alone) run for the next 500 yr into the future under the 
idealized condition that the current distribution of ice thick- 
ness prevails as is, thus imposing Ah(x, y, t) = 0 m for all 
t e [0, 500] yr. We find similar spatial patterns of bed uplift, 
as shown in Fig. 2a for the current uplift rate, in the future (cf. 
Supplement Fig. S7b); the total amount of GIA uplift at years 
AD 2100 and AD 2500 is in the respective ranges of about 
[—0.1, 0.5] (Fig. 4a) and [—0.5, 2.5] m (Fig. 4b). In such an 
idealized scenario of the unchanged future AIS, the notable 
features are that the peninsula, the whole of the WAIS, and 
coastal regions of the EAIS are likely to uplift, with the high- 
est uplift occurring around the Mt. Ellsworth territory, and 
that the interior of the EAIS may remain unaffected or sub- 
side minimally (cf. Supplement Fig. S7b). These are consis- 
tent with characteristic features of the employed ice loading 
history (cf. Supplement Fig. S2). 

Next, we compute model-average GIA solutions at years 
AD 2100 (Fig. 4c) and AD 2500 (Fig. 4d) by forcing the 
ISSM/GIA model by the future ice load changes alone. We 
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Figure 3. Estimates of the future Antarctic bed uplift. Model-average predictions for bed uplift at (a) AD 2100 and (b) AD 2500 under the 
R8 scenario. Associated standard deviations are shown in subplots (c) and (d), respectively. Calculations are made by forcing the calibrated 
ISSM/GIA model (Fig. 2) by the future ice loading (Supplement Fig. S3) obtained from the SeaRISE project (Bindschadler et al., 2013; 
Nowicki et al., 2013). See Supplement Fig. S6 for corresponding solutions for the bed uplift rate. Roughly similar spatial patterns are 
obtained, with the respective range of values [—5, 45] and [—15, 75] mmyr -1 at AD 2100 and AD 2500. 


impose A h(x,y,t) = 0m for all t e [—120, 0] kyr to ensure 
that the model is properly spun up. Alternative solutions are 
found by subtracting the GIA solutions associated with the 
past loading alone (Fig. 4a and b) from the corresponding 
final predictions (Fig. 3a and b). The resulting solutions are 
essentially identical to those shown in Fig. 4c and d, imply- 
ing that the principle of linear superposition holds (cf. Sup- 
plement Fig. S7). 

Comparing the estimates of GIA uplift associated with the 
past (Fig. 4a and b) and future ice loading alone (Fig. 4c and 
d) with those depicted in Fig. 3a and b, we find that the fu- 
ture ice loading dominates and that the contribution of long- 
term viscous creep associated with the past is only about one 
tenth the magnitude of the predicted GIA uplift. This sug- 
gests that the errors associated with the ice loading history 
may have minor consequences for the future predictions of 
the Antarctic bed topography, provided the differing scenar- 
ios properly sample the possible amplitude of the ice sheet 
loading. Significant changes in magnitude and timing of the 
loading history may, however, require that the different man- 
tle viscosity be employed in the ISSM/GIA model to predict 
the current uplift rates correctly (cf. Sect. 2.3). This, in turn, 
may yield the different contribution (not necessarily higher) 
of long-term viscous creep of the solid Earth to predictions 


of the future bed uplift. Nonetheless, it is important here to 
highlight the need for better constraining the DIHs, particu- 
larly in the recent past (during the past 1 kyr). 

3.3 Change in bed slope 

Strong spatial variability in the future GIA uplift (cf. 
Sect. 3.1) implies that the Antarctic bed slope will be mod- 
ulated in the future (e.g., Gomez et al., 2013; Konrad et al., 
2013). We compute current and future bedrock slopes (as- 
sociated with the model-average GIA solutions) following 

a h (x,y,t) = yjaljx, ~y, r) +a^(x, y, t), where a bx (x,y,t) 
and a\> y ( x,y,t ) are the v and y components of the bed slope, 
respectively. Figure 5b, for example, depicts the present-day 
bedrock slope of Antarctica. While this plot reveals the de- 
gree of bed steepness, it does not provide the information re- 
garding the aspect of slope. It is important to identify whether 
the bedrock has forward or reverse slope, particularly while 
evaluating the role of GIA uplift in the marine ice sheet in- 
stability (to be discussed later). We therefore plot the current 
bathymetry of the AIS in Fig. 5 a; in order to facilitate the dis- 
cussion, we only consider the areas with bedrock below sea 
level. Notice in the figure, for example, the blue color around 
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Figure 4. Role of past and future ice loading in the future bed topography. GIA uplift at (a) AD 2100 and (b) AD 2500 obtained by forcing 
the calibrated ISSM/GIA model (cf. Fig. 2.) by the past loading alone. Future DIHs are assumed to be zero. Corresponding solutions (model- 
average) at (c) AD 2100 and (d) AD 2500 associated with future ice loading alone. Past DIHs are assumed to be zero. For ease of comparison, 
we use the same color scale for each year. (Results in Figs. 4a and b are in the respective ranges of [—0.1, 0.5] and [—0.5, 2.5] m.) Also, see 
Supplement Fig. S7 where we replot some of these figures with a different color scale that better illustrates the spatial distribution. 


the interior of WAIS that illustrates the existence of reverse 
bed slope in those regions. 

We obtain the future bedrock topography by adding the 
current bedrock topography and the future GIA solution 
(Sect. 3.1). From the corresponding bedrock topography, 
we calculate the bed slope at AD 2100 and AD 2500. The 
changes in bedrock slope are then computed by subtracting 
the present-day slope (Fig. 5b) from the future bed slope. 
Sample results are shown in Fig. 5c for AD 2500. In the fig- 
ure, we generally notice the reduced bed slope (apart from 
a few localized regions with enhanced slope) around the 
Amundsen Sea Sector and the Amery Ice Shelf where large 
uplift is predicted. The reverse bed in these regions will thus 
have generally less steep slope in the future. Similar spatial 
patterns (but small magnitudes) of change in bed slope are 
obtained for AD 2100 as well (results not shown). 

Although the magnitudes of change are larger in the re- 
gions where the bed has experienced large uplift, such as 
in the Amundsen Sea Sector and in the Amery Ice Shelf, 
percent changes in bed slope are also significant around the 
Ronne and Ross ice shelves (Fig. 5d). The bedrock slopes 
beneath these ice shelves, for example, reduce by more than 
one percent at AD 2500. Under the circumstances when GL 
is to advance (this actually is the general case in the present 


context; cf. Sect. 4.2), change in bed slopes around the GL 
beneath the ice shelves may impact the magnitude of GL mi- 
gration (cf. Sect. 1.1 and Fig. la) and thus ice dynamics. 


4 Implications for ice sheet dynamics 

Our predictions for the future evolution of Antarctic bed 
topography may influence the future dynamics of AIS. In 
this section, we specifically quantify the potential effect of 
the predicted change in bed topography on the gravitational 
driving stress (Sect. 4.1), GL migration (Sect. 4.2), and ice 
surface velocities (Sect. 4.3). Because our model is not yet 
capable of computing IS/SE interactions with full dynamic 
feedbacks, it should be noted that some of the results pre- 
sented below are obtained by bootstrapping the relevant fu- 
ture bedrock topography. The general procedure includes the 
following. First, we consider the present-day settings (geo- 
metric setting and boundary conditions) of the AIS for our 
calculations. Next, we upgrade the geometry and relevant 
boundary conditions (e.g., basal friction while calculating ice 
surface velocities in Sect. 4.3) to account for the future GIA 
uplift and perform recalculations. Comparing corresponding 
results, we finally isolate the influence of GIA uplift. 


www.solid-earth.net/5/569/2014/ 


Solid Earth, 5, 569-584, 2014 


578 


S. Adhikari et al.: Antarctic ice dynamics and isostatically evolving bed 



Figure 5. Estimated change in Antarctic bed slope in the future. Present-day (a) bathymetry and (b) bed slope, a^(x, y ), of Antarctica. To 
facilitate discussions, only data in the range [—1000, 0] m are shown in Fig. 5a. (c-d) Model-average change in bed slope, Aa^ix, y), at 
AD 2500. Negative magnitudes of Aa^ix, y) imply that bedrock will have less steep slopes in the future. 


In order to minimize the potential compounding effects 
on the ice sheet dynamics, we keep the present-day settings 
(particularly, ice thickness and thermo-mechanical boundary 
conditions) fixed for the further calculations. We therefore 
advise caution in overinterpreting any individual results, ob- 
tained from the present-day settings of AIS perturbed by the 
predicted GIA uplift, for AD 2100 and 2500. For clarity, we 
revert to the general time convention (cf. Sect. 2.2) in order 
to stress that the following results illustrate the potential ef- 
fect of the predicted bed uplift, not the predictions of the total 
GIA effect, on the future dynamics of AIS. 

4.1 Gravitational driving stress 

Here, we discuss the GIA effects on the gravitational 
driving stress. We compute driving stress following 

T d (x,y, t ) = + xl y (x,y,t), where z^ix, y, t) 

and Td y (x,y,t) are the jc and y components of the driv- 
ing stress, respectively. For i = x,y, we define T& (x, y, t) = 
— p gh(x,y,t)a S i(x,y,t ) as the /th component of driving 
stress (e.g., Cuffey and Paterson, 2010), where g is the gravi- 
tational acceleration, a S i (= fy/z+a'b* ) is the i th component of 
surface slope, and 3 \h is the thickness gradient in the i th di- 
rection. In order to isolate the GIA effect, as noted earlier, we 
use the present-day ice thickness, i.e., h(x, y, t) = h(x, y, 0), 
and its gradients in all calculations. Hence, the change in bed 


slope due to the GIA uplift (Sect. 3.3) is responsible for mod- 
ulating the driving stress. 

Figure 6 shows the change in driving stress associated with 
the bed uplift predicted at years AD 2100 and AD 2500. In 
both cases (Fig. 6a and b), we notice small but similar trends 
of change in driving stress. Relatively large changes are pre- 
dicted at positions of larger bed uplift. Reduction in driving 
stress is evident around the Amery Ice Shelf, implying that 
the local surface slopes are likely to flatten in this region. 
Minor increments in driving stress are predicted in the area 
around Dome C. Complex patterns are predicted around the 
Amundsen Sea Sector; an extensive zone of reduced driving 
stress is surrounded by zones with enhanced driving stress 
(see particularly Fig. 6b). 

Theoretically, for the given distribution of ice thickness, 
changes in bed slope and surface slope (and, hence, the driv- 
ing stress) should generally be in phase for the topogra- 
phy with forward slope, but out of phase for cases with re- 
verse bed slope. Due to the complex nature of the Antarctic 
bathymetry, it is an arduously difficult task to find a robust 
and consistent relationship between changes in bed slope 
and driving stress, particularly around the Amundsen Sec- 
tor (compare, for example, Fig. 5c vs. Fig. 6b). Nonetheless, 
we generally find the reduction in local driving stress asso- 
ciated with the GIA uplift (compare, for example, Fig. 3b 
vs. Fig. 6b). Given the small order- of-magnitude predictions 
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Figure 6. Influence of GIA uplift on the gravitational driving stress. Change in driving stress, Ar&(x,y), due to the predicted bed uplift at 
(a) AD 2100 and (b) AD 2500. Calculations are made for the current distribution of ice thickness. Negative magnitudes imply that surface 
slopes flatten in the future. As the GIA solutions only perturb the ice-bedrock contact area, we mask out the ice shelves in the figures. 


for change in driving stress (i.e., several hundreds of Pascal 
only), it is important here to note that, as will also be shown 
in Sect. 4.3, minor changes in driving stress may affect the 
ice sheet dynamics substantially, as ice velocities are directly 
proportional to the driving stress by about a power of three 
(e.g., Cuffey and Paterson, 2010). 

4.2 GL migration 

In this section, we evaluate the effects of GIA uplift on the 
GL. We employ the simple hydrostatic equilibrium criterion 
to identify the transition points by seeking the floating ice 
thickness, hf(x, y, t ), such that 

hf(x,y,t) = — — b(x,y,t). (3) 

P 

Regions with ice thickness h(x, y, t) > hf(x, y, t) are as- 
sumed to be grounded and the rest to be floating. Because we 
use the present-day ice thickness, i.e., h(x, y, t) = h(x, y, 0), 
in all calculations, the extent of the current and future 
grounded ice (i.e., GL) is determined by the corresponding 
bathymetry, thus highlighting the influence of predicted GIA 
solutions. 

Using Eq. (3), we compute the extent of the grounded ice 
for t = 0, 100, and 500 yr. Changes in GL are then identi- 
fied by subtracting, in turn, the first solution from the latter 
two. Figure 7a, for example, shows the GL migration asso- 
ciated with the GIA uplift predicted for AD 2500. The mask 
shown in the figure primarily represents GL advance, imply- 
ing that more ice will be grounded in the fixture due to the 
GIA effects. However, we also predict the minor GL retreat 
in a few scattered locations, particularly along the Wilkes 
Land coast, where the bedrock is generally predicted to sub- 
side partly due to the large snow accumulation simulated 
under the chosen climate change scenario (Nowicki et al., 
2013). Figure 7b, for example, depicts the mask of GL retreat 


around the Shackleton Ice Shelf. Note that we also obtain 
similar but less extensive migration in GL associated with 
the GIA uplift predicted for AD 2100 (results not shown); 
there is however no evidence of GL retreat in this case. 

Based on the measured ice surface velocities (Rignot et al., 
2011), we locate more than 2700 ice flowlines (cf. Supple- 
ment Fig. S8) to quantify the magnitude of GL migration 
(mostly advance) associated with the predictions of GIA up- 
lift at AD 2500. Results are shown for two important regions, 
namely the Ronne (Fig. 7c) and Ross (Fig. 7d) ice shelves. 
Figures reveal that the GL may advance by more than 1 00 km 
in these two ice shelves. Significant GL advance (by tens of 
km) is also predicted in the Amery Ice Shelf, Amundsen Sec- 
tor, Larsen Ice Shelf, Brunt Ice Shelf, and in other regions. 
In a few locations, e.g., the Shackleton Ice Shelf (Fig. 7b), as 
noted earlier, we predict a minor retreat (by <10 km) of the 
GF. 

Although the primary control of GF migration in our cal- 
culations is bedrock elevation (Eq. 3), it should be noted that 
the bedrock slope plays an equally important role (e.g., Lin- 
gle and Clark, 1985; Gomez et al., 2010), as summarized in 
Sect. 1.1. Extensive advancement in GL associated with the 
GIA uplift is generally consistent with what is expected when 
slope reduces in the reverse bedslope topography (compare 
Fig. 5a vs. c, for example, around the Ronne Ice Shelf). 
For bedrock with forward slopes, however, the advancement 
in GL can be explained by the GIA-induced increment in 
bed slopes (compare Fig. 5a vs. c, for example, around the 
Getz Ice Shelf). The minor GF retreat predicted, for exam- 
ple, in the Shackleton Ice Shelf is associated with the flat- 
tening (Fig. 5c) of the forward bed slope (Fig. 5a). Although 
we are able to show a systematic relationship between the 
bathymetry, change in bed slope, and the direction of the 
GF migration in a few cases, it is complicated to provide 
such one-to-one relationships consistently over the entire ice 
sheet. Nonetheless, the results indicate that the GIA effects 
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Figure 7. Influence of GIA uplift on the GL. (a) The mask of GL migration associated with the GIA solution at AD 2500. Calculations are 
based on the hydrostatic equilibrium criterion for the current distribution of ice thickness. Cyan depicts the extent of present-day grounded 
ice. Red shows the GL migration (mostly advance) due to the GIA uplift. Minor retreats in GL are predicted in a few areas along Wilkes 
Land in the EAIS; their extents are limited to one or two mesh elements, i.e., < 10 km. (b) Mask of example retreat around the Shackleton Ice 
Shelf, as seen by zooming into the small blue box in Fig. 7a. Other large boxes in Fig. 7a enclose two important regions that are magnified: 
(c) the Ronne Ice Shelf and the Bellingshausen Sea Sector; and (d) the Ross Ice Shelf. Color codes illustrate the magnitude of GL advance 
measured along the ice flowlines. 


generally support the extension of grounded ice (i.e., GL ad- 
vance) in the future, thereby promoting the stability in the 
marine portions of the ice sheet that rest on a reverse bed 
slope. 

4.3 Ice surface velocities 

Finally, we analyze the influence of GIA uplift on the ice sur- 
face velocities. By solving the quasi-static thermomechanical 
problem of ice flow, we calculate the englacial velocity field 
of the AIS with and without GIA effects. We assume that 
higher-order physics based on the equations governing mass 
and momentum conservation (Pattyn, 2003) together with the 
constitutive relations for isotropic ice (Glen, 1955) describe 
the internal creep deformation of ice, and that a viscous law 
of friction (e.g., MacAyeal, 1993) governs basal sliding. We 
rely on a steady-state thermal problem identical for all sim- 
ulations, based on present-day conditions. For simplicity, we 
do not consider the possibility of till deformation underneath 
the ice sheet. A description of ice rheology (Glen, 1955) and 
other common assumptions related to ice flow modeling can 
be found elsewhere (e.g., Cuffey and Paterson, 2010). 


For the present-day setting of the AIS, we solve this prob- 
lem of ice dynamics through diagnostic simulation of the 3-D 
ice flow capability of ISSM, satisfying a number of boundary 
conditions (e.g., Larour et al., 2012b). We impose (i) a stress- 
free condition at the ice-atmosphere interface, (ii) water 
pressure directing towards the ice sheet at the ice-ocean (pe- 
ripheral) interface, and (iii) a sliding condition governed by 
the basal friction at the ice-bed interface. Zero friction is ap- 
plied at the base of the ice shelf (i.e., free-floating condition), 
while basal friction under the grounded ice is inferred from 
InSAR (Interferometric Synthetic Aperture Radar) based sur- 
face velocities (Rignot et al., 2011) using a data assimilation 
technique (e.g., Morlighem et al., 2010). The basal friction 
pattern is similar to the one described in Morlighem et al. 
(2013). 

We re-run the simulation for two additional cases, asso- 
ciated with the predictions of GIA uplift at years AD 2100 
and AD 2500. In each case, we upgrade the bedrock and 
surface elevations of the grounded ice; the extent of the 
grounded ice (GLs) is also updated (cf. Sect. 4.2). For float- 
ing ice as well, we upgrade both bed and surface elevations 
so that all floating ice is in hydrostatic equilibrium, which 
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Figure 8. Influence of GIA uplift on the ice surface velocities. Change in ice surface velocity, A u(x, y), due to the predicted bed uplift at 
(a) AD 2100 and (b) AD 2500. Calculations are made by running the diagnostic simulation of the 3-D ice flow (higher-order) capability 
of ISSM. Other model setup and boundary conditions are consistent with those of the SeaRISE control experiment (cf. Bindschadler et al., 
2013). A systematic reduction in velocity, which can be attributed partly to the reduction in z^(x, y) (around the Amundsen Sector and the 
Amery Ice Shelf; cf. Fig. 6) and partly to the GL advance (particularly in the large ice shelves; cf. Fig. 7), indicates the stabilizing effects of 
GIA uplift on the future dynamics of AIS. 


ensures continuity of the driving stress at the GL. Apply- 
ing the same boundary conditions discussed above (except 
in the newly grounded or floating areas), we compute the 
englacial velocity for both cases associated with the future 
GIA uplift. In areas presently floating that become grounded 
at t = 100 and 500 yr (cf. Fig. 7a), we update basal friction 
assuming that it is roughly equal to the gravitational driving 
stress (Morlighem et al., 2013). Similarly, we impose a free- 
floating condition in areas presently grounded that float in the 
future (cf. Fig. 7b). 

For a given vertical profile of an ice sheet, the maximum 
velocity is always observed at the surface. We therefore place 
our emphasis upon the GIA influence on the ice surface ve- 
locities. Using the simulation results discussed above, we 
compute the GIA-induced change in surface velocity asso- 
ciated with the predictions of GIA uplift at years AD 2100 
(Fig. 8a) and AD 2500 (Fig. 8b). In both cases, we find 
similar patterns of change in ice surface velocity. Although 
the predicted changes are small, about two to three orders 
of magnitude smaller than the surface velocities themselves 
(Rignot et al., 2011), a systematic reduction in velocity is 
evident around the sheet-shelf margins. This suggests that 
the GIA effects generally contribute to decelerating the flow 
speed across the GL, and hence promote stability in the ma- 
rine portions of the ice sheet. Note that we mask out the ice 
shelves in our figures, for we have no intention of making 
predictions in the ice shelves. 

The predicted changes in surface velocity for the grounded 
ice can be interpreted as the combined effects of changes 
in driving stress (cf. Sect. 4.1 and Fig. 6) and the GL (cf. 
Sect. 4.2 and Fig. 7) associated with the GIA uplift. Around 
the Ross and Ronne ice shelves, the GIA-induced reduction 
in surface velocity is consistent with the GL advance. In other 


regions, e.g., the Amundsen Sea Sector and the Amery Ice 
Shelf, predicted reduction in ice velocity can be attributed 
partly to the GL advance and partly to the reduced driving 
stress. All in all, the effects of GIA on several aspects of ice 
dynamics (e.g., driving stress, GL, and ice surface velocities) 
are consistent in that the GIA promotes systematic stability 
in marine portions of the AIS in the future. 


5 Conclusions 

This study has examined the interplay between the ice sheet 
evolution and the solid Earth responses for the AIS. First, 
we compute the future uplift of the Antarctic bedrock using 
the calibrated ISSM/GIA model driven by the inferred and 
predicted evolution of the ice sheet. Next, we evaluate how 
such a response of the solid Earth impacts AIS dynamics. 

Our calculations are based on several approximations of 
model physics and numerics; it is important to highlight 
some of these here. The GIA model describes a simple two- 
layer representation of the solid Earth; the model and ma- 
terial viscoelastic parameters are kept constant spatially and 
the viscosity is Newtonian. Our ice sheet model solves the 
quasi-static thermomechanical flow problem for higher-order 
mechanics; the GLs are determined by the hydrostatic equi- 
librium criterion. A more comprehensive exploration of the 
positive or negative IS/SE feedbacks is warranted in the fu- 
ture. There is much to be learned from GIA models that em- 
ploy additional GPS data, possibly right in the heart of the 
Amundsen Sea Sector, where viscoelastic uplift rates may 
approach 40mmyr -1 (Groh et al., 2012). Our computations 
of ice loading after the present-day rely on several ice sheet 
models driven by the melt-dominated forcing under the R8 
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scenario. The computed model-average GIA response pro- 
vides our assessment of its impact on the ice sheet dynam- 
ics. While there are limitations to the data and methods em- 
ployed, our research brings us to the following two important 
conclusions of broader interest. 

First, the short-term viscoelastic response of the solid 
Earth to the future ice load change, rather than its long-term 
viscous response to the past loading, governs the future evo- 
lution of the Antarctic bed topography. The magnitude and 
spatial variability in the future bed uplift are therefore de- 
termined by the nature of the future evolution of the AIS. 
A larger uplift is expected, for example, where the ice sheet 
loses more mass, while its far-held consequences seem to 
involve a relatively small amplitude subsidence. Our calcu- 
lations suggest that the Antarctic bed may rise by a few me- 
ters and a few tens of meters around the WAIS, particularly 
the Amundsen Sea Sector and the Amery Ice Shelf, at years 
AD 2100 and AD 2500, respectively. Minor subsidences of 
about one meter and a few meters are predicted along the 
Wilkes Land coast at the respective times, partly caused by 
the net accumulation in the climate scenario runs (Nowicki 
et al., 2013). The interior of the EAIS is likely to remain un- 
changed. 

Second, a pervasive and large uplift predicted in the in- 
terior of the WAIS, a substantially marine-based ice sheet, 
has particular significance because it generally corresponds 
to the battening of the reverse bed slope. This drives the GL 
forward and consequently promotes the stability of the ice 
sheet. Our calculations, based on the present-day setting of 
the AIS perturbed by the predicted GIA uplift, reveal that the 
GL may advance by more than 1 00 km in the Ross and Ronne 
ice shelves due to the predicted GIA uplift for AD 2500. This 
may reduce the future ice surface velocities across the GLs 
by several tens of meter per annum. 

The conclusions summarized above indicate a negative 
feedback between the ice sheet evolution and the solid Earth 
response for the marine ice sheet. For areas with reverse bed 
slope, for example, loss in ice mass battens the bed and drives 
the GL forward and hence decelerates the rate of mass loss. 
Although our model is capable of illustrating this mechanism 
systematically, accurate quantibcation of its signibcance re- 
quires a dynamically coupled IS/SE model including the in- 
tricate details of ice sheet stability to breakup. This negative 
feedback is consistent with the ice sheet and sea level simu- 
lations computed by Gomez et al. (2010, 2013) wherein loss 
in ice mass reduces the local sea level due to self-gravitation 
and hence decelerates the rate of mass loss. For accurate sim- 
ulations of the AIS on centennial timescales under the rea- 
sonable climate change scenarios, both the solid Earth and 
sea level changes proximal to the GL may need to be prop- 
erly accounted for. 


The Supplement related to this article is available online 
at doi:10.5194/se-5-569-2014-supplement. 
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